SEA rsfcMRI - Longitudinal Descriptives
John Flournoy
2020-10-14
fslong <- FALSE ### CHANGE FILE NAME TOO!
library(ggplot2)
library(data.table)
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(tidyr)
library(showtext)## Loading required package: sysfonts## Loading required package: showtextdb#{.tabset}
font_add_google("Didact Gothic", "Didact Gothic")
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)
#https://www.instagram.com/p/CC3EvmLgt2m/?utm_source=ig_web_copy_link
apal <- paste0('#',c('2C2B2B', 'F9773B', 'FFEA8C', '1389E6', 'D2E5E7'))
jftheme <- theme_minimal() +
theme(text = element_text(family = 'Didact Gothic', size = 14),
panel.background = element_rect(fill = apal[[5]], size = 0, color = apal[[2]]),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
strip.background = element_rect(fill = apal[[3]], size = 0),
strip.text = element_text(color = '#222222'),
axis.text = element_text(color = apal[[1]]), axis.title = element_text(color = apal[[1]]))
knitr::read_chunk('collect_data.R')
knitr::read_chunk('fit_models.R')if(fslong){
cat('\n\n# FreeSurfer Longitudinal Registration\n\n')
} else {
cat('\n\n# fMRI-prep Basic Registration\n\n')
}fMRI-prep Basic Registration
Load the data
Starting with the Schaefer 400-parcel, 7 network atlas.
Retain just within-network connectivity. Also, Fisher-z transform the correlations for analyses. One small detail that is important here is that we keep network connectivity for same-hemisphere parcels only.
library(data.table)
library(dplyr)
library(tidyr)
#set variables for `collect_data`
if(fslong){
fname <- 'sea_rsfc_fslong_schaefer400x7.RDS'
data_dir <- '/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/SEA_REST/derivatives/fmriprep-20.1.1/xcpengine-default/'
subs <- sprintf('sub-10%02d', c(1:16,18:31))
sess <- sprintf('%02d', 1:10)
subses <- expand.grid(sess, subs)
files <- paste0(data_dir, '/', subses[,2], subses[,1], '/fcon/schaefer400x7/', subses[,2], subses[,1], '_schaefer400x7.net')
file_id <- paste0(subses[,2], '_', subses[,1])
} else {
fname <- 'sea_rsfc_schaefer400x7.RDS'
data_dir <- '/net/holynfs01/srv/export/mclaughlin/share_root/stressdevlab/SEA_BIDS/derivatives/1.4.1-final/xcpengine-default'
subs <- sprintf('sub-10%02d', c(1:16,18:31))
sess <- sprintf('ses-%02d', 1:10)
subses <- expand.grid(sess, subs)
files <- paste0(data_dir, '/', subses[,2], '/', subses[,1], '/fcon/schaefer400x7/', subses[,2], '_', subses[,1], '_schaefer400x7.net')
file_id <- paste0(subses[,2], '_', subses[,1])
}
library(data.table)
library(tidyr)
if(is.na(ncores <- as.numeric(Sys.getenv('SLURM_CPUS_ON_NODE')))){
ncores <- 3
NOTSLURM <- TRUE
message('No environment variable specifying number of cores. Setting to ', ncores, '.')
} else {
NOTSLURM <- FALSE
message('Number of cores set to ', ncores)
}
setDTthreads(ncores)
if(file.exists(fname)){
adt_labels <- readRDS(fname)
schaefer_400_7_net_ids <- readRDS('schaefer_ids.RDS')
} else {
message('Creating file list...')
files_list <- lapply(files, function(x) ifelse(file.exists(x), x, ''))
message('Reading rsfc files...')
if(ncores > 1){
message('Using ', ncores, ' cores to read files.')
cl <- parallel::makePSOCKcluster(max(c(ncores - 1), 1))
adt_list <- unlist(parallel::parLapply(cl = cl, split(files, 1:ncores), function(part){
lapply(part, function(f){
if(file.exists(f)){
data.table::fread(f, skip = 2, col.names = c('node1', 'node2', 'r'))
} else {
data.table::data.table()
}
})
}), recursive = FALSE)
try(parallel::stopCluster(cl))
} else {
adt_list <- lapply(files, function(f){
if(file.exists(f)){
data.table::fread(f, skip = 2, col.names = c('node1', 'node2', 'r'))
} else {
data.table::data.table()
}
})
}
names(adt_list) <- file_id
message('Combining data, assigning labels...')
adt <- rbindlist(adt_list, idcol = 'id')
#https://github.com/ThomasYeoLab/CBIG/tree/master/stable_projects/brain_parcellation/Schaefer2018_LocalGlobal/Parcellations/MNI
schaefer_400_7_net_ids <- fread('Schaefer2018_400Parcels_7Networks_order.txt',
col.names = c('node1', 'network1', 'V1', 'V2', 'V3', 'V4'))
schaefer_400_7_net_ids[, c('V1', 'V2', 'V3', 'V4') := NULL]
schaefer_400_7_net_ids <- schaefer_400_7_net_ids %>%
tidyr::extract(network1, c('network1', 'roi1'), '7Networks_([RL]H_.*)_(\\d+)')
setkey(adt, node1)
setkey(schaefer_400_7_net_ids, node1)
adt_labels1 <- schaefer_400_7_net_ids[adt]
setnames(schaefer_400_7_net_ids, c('node2', 'network2', 'roi2'))
setkey(adt_labels1, node2)
setkey(schaefer_400_7_net_ids, node2)
adt_labels <- schaefer_400_7_net_ids[adt_labels1]
adt_labels[, c('id', 'sess') := tstrsplit(id, '_', fixed = TRUE)]
message('Saving RDS file to: ', fname)
saveRDS(adt_labels, fname)
message('Saving Schaefer ID data table to schaefer_ids.rds')
saveRDS(schaefer_400_7_net_ids, 'schaefer_ids.RDS')
}
sea_schaefer400_7_withinnet <- copy(adt_labels)
sub_regex <- '[LR]H_[A-Za-z]+_*(.*)'
net_regex <- '[LR]H_([A-Za-z]+)_*(?:.*)'
hem_regex <- '([LR]H)_[A-Za-z]+_*(?:.*)'
sea_schaefer400_7_withinnet_nets <- distinct(sea_schaefer400_7_withinnet, network1)
sea_schaefer400_7_withinnet_nets[, sub1 := gsub(sub_regex, '\\1', network1)]
sea_schaefer400_7_withinnet_nets[, net1 := gsub(net_regex, '\\1', network1)]
sea_schaefer400_7_withinnet_nets[, hem1 := gsub(hem_regex, '\\1', network1)]
setkey(sea_schaefer400_7_withinnet, network1)
setkey(sea_schaefer400_7_withinnet_nets, network1)
sea_schaefer400_7_withinnet <- sea_schaefer400_7_withinnet_nets[sea_schaefer400_7_withinnet]
setnames(sea_schaefer400_7_withinnet_nets, c('network2', 'sub2', 'net2', 'hem2'))
setkey(sea_schaefer400_7_withinnet, network2)
setkey(sea_schaefer400_7_withinnet_nets, network2)
sea_schaefer400_7_withinnet <- sea_schaefer400_7_withinnet_nets[sea_schaefer400_7_withinnet]
sea_schaefer400_7_withinnet <- sea_schaefer400_7_withinnet[net1 == net2 & hem1 == hem2]
sea_schaefer400_7_withinnet[, c('network1', 'network2') := NULL]
sea_schaefer400_7_withinnet[, z := atanh(r)]
sea_schaefer400_7_withinnet[, H := dplyr::case_when(hem1 == 'RH' ~ -1,
hem1 == 'LH' ~ 1)]
if(!dim(distinct(sea_schaefer400_7_withinnet, net1, net2, hem1, hem2))[[1]] == 14){
stop('Incorrect number of networks')
}
if(FALSE){
adt_labels_old <- readRDS('sea_rsfc_schaefer400x7.RDS')
missing_files_csv <- data.table(file = files, does_not_exist = files_list == '')
View(missing_files_csv[does_not_exist == TRUE])
readr::write_csv(data.table(file = files, does_not_exist = files_list == ''),
'~/sea_rsfc-missing_fcon_output.csv')
a <- data.table(file = files, does_not_exist = files_list == '')
a[, fcon_missing := lapply(gsub('(.*/fcon)/.*', '\\1', file), function(x) !file.exists(x))]
a[, missmatch := does_not_exist != fcon_missing ]
a[(missmatch), file]
}Check No. of RS Features
Using the number nodes are in each network allows computation of the expected number of rsfc features, which allows a comparison to the observed number of features to find instances where certain connections are missing.
#Check to make sure we're getting the number of connectivity features per network we expect
setkey(sea_schaefer400_7_withinnet_nets, 'network2')
setkey(schaefer_400_7_net_ids, 'network2')
nodes_per_net <- sea_schaefer400_7_withinnet_nets[schaefer_400_7_net_ids][, .N, by = c('net2', 'hem2')]
nodes_per_net[, expected_Nfeatures := N * (N - 1) / 2]
connectivity_feature_check <- sea_schaefer400_7_withinnet[, .(Nfeatures = .N), by = c('net2', 'hem2', 'id', 'sess')][nodes_per_net[, !'N'], on = c('net2', 'hem2')]
setnames(connectivity_feature_check, c('net2', 'hem2'), c('nework', 'hemisphere'))
dt_options <- list(rownames = FALSE,
filter = 'top',
class = 'cell-border stripe',
extensions = 'Buttons',
options = list(dom = 'Bfrtip', buttons = c('csv')))
do.call(DT::datatable,
c(list(connectivity_feature_check[Nfeatures != expected_Nfeatures],
caption = 'Mismatch of observed and expected number of rsfc features'),
dt_options))do.call(DT::datatable,
c(list(connectivity_feature_check[Nfeatures != expected_Nfeatures][, .N, by = c('id', 'sess')],
caption = 'Number of mismatches for subject sessions'),
dt_options))Node coverage
Each node in a network should be connected to N - 1 other nodes in the network, where N is the number of total nodes in the network.
Mean connections-per-node
node_check <- sea_schaefer400_7_withinnet[, c('node1', 'node2', 'id', 'sess', 'net2', 'hem2')]
node_check[, idx := 1:.N, by = c('id', 'sess', 'net2', 'hem2')]
node_check[, c('maxnode', 'minnode') := .(max(node1, node2), min(node1, node2)), by = c('node1', 'node2')]
node_check_l <- melt(node_check, measure.vars = c('maxnode', 'minnode'), value.name = 'node')
node_info <- node_check_l[, .N,
by = c('id', 'sess', 'net2', 'node', 'hem2')][, list(total = max(N)),
by = c('node', 'net2', 'hem2')]
net_info <- node_info[, .(total = unique(total)), by = c('net2', 'hem2')]
node_check_subj <- node_check_l[, .N, by = c('sess', 'id', 'node', 'net2', 'hem2')]
node_check_subj <- net_info[node_check_subj, on = c('net2', 'hem2')]
node_check_subj_summary <- node_check_subj[, .(observed = mean(N),
expected = mean(total),
proportion = sum(N) / sum(total)),
by = c('node', 'net2', 'hem2')]
node_check_subj_summary[, network := paste(net2, hem2, sep = '_')]
node_check_subj_summary_l <- melt(node_check_subj_summary, measure.vars = c('observed', 'expected'))
node_check_subj_summary_l[, variable := factor(variable, levels = c('expected', 'observed'))]
ggplot(node_check_subj_summary_l,
aes(x = node, y = value)) +
geom_line(aes(fill = variable, color = variable), position = position_identity()) +
scale_color_manual(values = apal[c(1,2)], aesthetics = c('color', 'fill')) +
facet_wrap(~ network, ncol = 2, scales = 'free') +
jftheme +
labs(y = 'Mean number of connections from node', x = 'Node')
Proportion connections-per-node
ggplot(node_check_subj_summary,
aes(x = node, y = proportion)) +
geom_line(color = apal[[4]]) +
facet_wrap(~ network, ncol = 2, scales = 'free') +
jftheme +
labs(y = 'Proportion of expected connections from node', x = 'Node')
Mean across sessions by subject
node_check_sess_summary <- node_check_subj[, .(mean = mean(N),
proportion = sum(N) / sum(total)),
by = c('id', 'net2', 'hem2', 'sess')]
node_check_sess_summary[, sess_num := as.numeric(gsub('ses-(\\d+)', '\\1', sess))]
ggplot(node_check_sess_summary,
aes(x = sess_num, y = mean)) +
geom_line(aes(group = id, color = id), alpha = .5) +
facet_grid(net2 ~ hem2, scales = 'free') +
jftheme +
labs(y = 'Mean variability per node', x = 'Session')
Proportion across sessions by subject
node_check_sess_summary <- node_check_subj[, .(mean = mean(N),
proportion = sum(N) / sum(total)),
by = c('id', 'net2', 'hem2', 'sess')]
node_check_sess_summary[, sess_num := as.numeric(gsub('ses-(\\d+)', '\\1', sess))]
ggplot(node_check_sess_summary,
aes(x = sess_num, y = proportion)) +
geom_line(aes(group = id, color = id), alpha = .5) +
facet_grid(net2 ~ hem2) +
jftheme +
labs(y = 'Proportion variability per node', x = 'Session')
Another proportion across subjects
node_check_sess_summary[, proportion_scaled := proportion / .N * 10, by = c('hem2', 'net2', 'id')]
id_props <- node_check_sess_summary[, list(prop_mean = mean(proportion, na.rm = TRUE)), by = c('id')]
node_check_sess_summary[, id_sorted := factor(id, levels = id_props$id[order(id_props$prop_mean)])]
ggplot(node_check_sess_summary, aes(x = id_sorted, y = proportion_scaled)) +
geom_col(aes(group = sess, fill = proportion, color = proportion),
position = 'stack') +
scale_fill_gradientn(colors = apal[c(2,1)], values = c(0, 1),
aesthetics = c('fill', 'color')) +
facet_grid(net2 ~ hem2) +
jftheme +
theme(axis.text.x = element_text(angle = 360-65, hjust = 0, size = 6))
Overall proportions table
setorder(id_props, 'prop_mean')
id_props[, `Mean prop.` := sprintf('%0.2f', prop_mean)]
do.call(DT::datatable,
c(list(id_props[, c('id', 'Mean prop.')],
caption = 'Average proportion coverage'),
dt_options))Some QC
Density of functional connectivity (pearson correlation) for each participant, for each subject, overlayed on the density for all sessions and participants.
#Generate group-level density
agg_mean <- mean(sea_schaefer400_7_withinnet$z)
agg_sd <- sd(sea_schaefer400_7_withinnet$z)
scalez <- (sea_schaefer400_7_withinnet$z - agg_mean) / agg_sd
density_agg <- density(scalez)
sea_schaefer400_densplot_agg <- data.table(x = density_agg$x,
x_on_z = density_agg$x * agg_sd + agg_mean,
y = density_agg$y)
#Generate per-participant-session densities
sea_schaefer400_MSD <- sea_schaefer400_7_withinnet[, list(mean = mean(z),
sd = sd(z)),
by = c('id', 'sess')]
sea_schaefer400_densplot <- sea_schaefer400_MSD[sea_schaefer400_7_withinnet,
on = c('id', 'sess')]
sea_schaefer400_densplot[, scalez := (z - mean)/sd]
sea_schaefer400_densplot <-
sea_schaefer400_MSD[sea_schaefer400_densplot[, list(x = density(scalez)$x,
y = density(scalez)$y),
by = c('id', 'sess')],
on = c('id', 'sess')][, x_on_z := sd*x + mean]u_ids <- unique(sea_schaefer400_densplot$id)
for(an_id in u_ids){
cat(paste0('\n\n## ', an_id, '\n\n'))
hplot <- ggplot(sea_schaefer400_densplot_agg, aes(x = tanh(x_on_z), y = y)) +
geom_ribbon(ymin = 0, aes(ymax = y, fill = 'Group', color = 'Group'),
alpha = .5) +
geom_ribbon(data = sea_schaefer400_densplot[id == an_id, ],
ymin = 0, aes(ymax = y, fill = 'ID', color = 'ID'),
alpha = .5) +
scale_fill_manual(aesthetics = c('color','fill'), breaks = c('Group', 'ID'), labels = c('Group', 'Participant'), values = apal[c(2,4)], name = 'Data') +
facet_wrap(~sess, ncol = 2) +
labs(x = 'correlation', y = 'density') +
coord_cartesian(xlim = c(-1, 1), ylim = c(0, .5)) +
jftheme
print(hplot)
}sub-1001
sub-1002
sub-1003
sub-1004
sub-1005
sub-1006
sub-1007
sub-1008
sub-1009
sub-1010
sub-1011
sub-1012
sub-1013
sub-1014
sub-1015
sub-1016
sub-1018
sub-1019
sub-1020
sub-1021
sub-1022
sub-1023
sub-1024
sub-1025
sub-1026
sub-1027
sub-1028
sub-1029
sub-1030
sub-1031
Modeling
We are interested in stability and variability for each network. For each participant, we have multiple sessions, and within each session, we have multiple parcel-parcel pairs that provide information about the within-network connectivity. We can use the fact that we have multiple indicators of within-network connectivity to estimate the between-person variability, as well as the within-person (session-to-session) variability as distinct from measurement error (we assume that each parcel-parcel functional connectivity estimate in a network is an estimate of that network's connectivity)*.
* This assumption means that we essentially treat differences in the FC between one pair and another as measurement error.
We can compute an ICC that describes both within-person and between-person variability by using a 3 level model. First, I subset the data for a single network. I then build an intercept-only model, allowing the intercept to vary by participant-ID, and by session within each ID. The intercept is effectively the mean across all intrahemispheric parcel-pairs.
\[ \begin{align} z &= \beta_{0jk} + \epsilon_{ijk} \\ \beta_{0jk} &= \pi_{00k} + \nu_{0jk} \\ \pi_{00k} &= \gamma_{000} + \upsilon_{00k} \end{align} \]
\[ \begin{align} z = \gamma_{000} + \upsilon_{00k} + \nu_{0jk} + \epsilon_{ijk} \end{align} \] Where \(z\) is the measure of functional connectivity, and \(i\) indexes observations within each session, \(j\), for each participant \(k\). This means that we are able to estimate variance in per-person deviations (\(\upsilon_{00k}\)) from the network mean-connectivity (\(\gamma_{000}\)), per-session-deviations (\(\nu_{0jk}\)) from those per-person deviations, and the residual not explained by variability over persons or person-sessions (\(\epsilon_{ijk}\)).
Testing these out
I want to make sure that these models are correct, so I'll compare the variance portions, and total variance, between the two models.
library(lme4)
networks <- unique(sea_schaefer400_7_withinnet$net1)
this_net <- networks[[2]]
this_net_dt <- sea_schaefer400_7_withinnet[net1 == this_net]
fslong_insert <- ifelse(fslong, '_fslong', '')
model_dir <- 'models'
test_model_list <- list(test_2l = list(fn = file.path(model_dir, paste0('test_noh', fslong_insert, '_2l.RDS')),
form = z ~ 1 + (1 | id)),
test_3l = list(fn = file.path(model_dir, paste0('test_noh', fslong_insert, '_3l.RDS')),
form = z ~ 1 + (1 | id/sess)))
cl <- parallel::makePSOCKcluster(max(c(ncores - 1), 1))
test_model_fits <- parallel::parLapply(cl = cl, test_model_list, function(model_list, d){
library(lme4)
if(!file.exists(model_list[['fn']])){
f <- model_list[['form']]
fit <- lmer(f, data = d)
saveRDS(fit, model_list[['fn']])
} else {
fit <- readRDS(model_list[['fn']])
}
return(fit)
}, d = this_net_dt)
try(parallel::stopCluster(cl))three_level <- test_model_fits$test_3l
summary(three_level)## Linear mixed model fit by REML ['lmerMod']
## Formula: z ~ 1 + (1 | id/sess)
## Data: d
##
## REML criterion at convergence: 186542
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.8941 -0.6869 -0.0774 0.6111 6.0356
##
## Random effects:
## Groups Name Variance Std.Dev.
## sess:id (Intercept) 3.194e-03 0.056516
## id (Intercept) 5.240e-06 0.002289
## Residual 8.027e-02 0.283317
## Number of obs: 587200, groups: sess:id, 289; id, 30
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 0.219506 0.003371 65.11lll_varcorr <- VarCorr(three_level)
report_df <- data.frame(
stat = c('id_var', 'sess_var', 'resid_var', 'total_var'),
three_level = c(lll_varcorr$id.1['(Intercept)','(Intercept)'],
lll_varcorr$sess.id.1['(Intercept)','(Intercept)'],
sigma(three_level)^2,
sum(unlist(lapply(lll_varcorr, diag))) + sigma(three_level)^2))
report_df$three_level_perc <-
sprintf('%.1f%%', report_df$three_level/report_df$three_level[[4]]*100)
report_df$three_level_RE_perc <-
c(sprintf('%.1f%%', report_df$three_level[1:2]/sum(report_df$three_level[1:2])*100), '-', '-')
knitr::kable(report_df, digits = 5)| stat | three_level | three_level_perc | three_level_RE_perc |
|---|---|---|---|
| id_var | 0.08027 | 96.2% | 49.0% |
| sess_var | 0.08347 | 100.0% | 51.0% |
| resid_var | 0.08027 | 96.2% | - |
| total_var | 0.08347 | 100.0% | - |
Fit for each network
networks <- unique(sea_schaefer400_7_withinnet$net1)
netlist <- lapply(networks, function(this_net){
return(list(network = this_net, d = sea_schaefer400_7_withinnet[net1 == this_net]))
})
if(NOTSLURM){
message('Not running as batch job. Will only attempt to read saved model fits.')
}
cl <- parallel::makePSOCKcluster(max(c(ncores - 1, 1)))
model_fits <- parallel::parLapply(cl = cl, netlist, function(anet, fslong, NOTSLURM){
this_net <- anet[['network']]
this_net_dt <- anet[['d']]
model_dir <- 'models'
fslong_name <- ''
if(fslong)
fslong_name <- 'fslong-'
model_fit_fn <- file.path(model_dir, paste0('lmer_fit-3l-', fslong_name, this_net, '.RDS'))
library(lme4)
if(!file.exists(model_fit_fn) & NOTSLURM){
stop('Model fit has not been saved. Will not estimate. Stopping.')
}
if(!file.exists(model_fit_fn)){
message('Fitting model for ', this_net)
fit <- lmer(z ~ 1 + (1 | id/sess), data = this_net_dt)
saveRDS(fit, model_fit_fn)
} else {
fit <- readRDS(model_fit_fn)
}
return(fit)
}, fslong = fslong, NOTSLURM = NOTSLURM)
try(parallel::stopCluster(cl))proportion_variance_tables <- lapply(model_fits, function(model_fit){
vc <- VarCorr(model_fit)
id_v <- vc$id['(Intercept)','(Intercept)']
idsess_v <- vc$`sess:id`['(Intercept)','(Intercept)']
s_2 <- sigma(model_fit)^2
total_RE <- sum(unlist(lapply(vc, diag)))
other_RE <- total_RE - id_v - idsess_v
total <- total_RE + s_2
rez <- data.frame(source = rep(c('ID', 'ID/Session', 'Other RE', 'Residual'),2),
out_of = factor(c(rep('total', 4), rep('RE', 4)), levels = c('total', 'RE')),
percent = c(c(id_v, idsess_v, other_RE, s_2)*100 / total,
c(id_v, idsess_v, other_RE, NA)*100 / total_RE))
if(other_RE == 0){
rez <- rez[rez$source != 'Other RE',]
}
return(rez)
})The plots on the left, below, show model-expected functional connectivity pearson correlations for each participant, for each session, overlayed on the raw data (one point per parcel, per session, per participant). A line for each participant's model-expected mean across all sessions is overlayed on this. The right plot shows proportions of variance accounted for by the random effect estimates as a proportion of both total variance, and total random effects (RE) variance. The total RE variance includes individual means (ID), individual means per session (ID/Session), and the difference in connectivity between right and left hemispheres (we treat this as a nuisance variable, essentially).
library(patchwork)
plot_percents <- function(percent_table){
ggplot(percent_table, aes(y = percent, x = out_of, fill = source)) +
geom_col(position = 'stack') +
scale_fill_manual(breaks = c('ID', 'ID/Session', 'Other RE', 'Residual'),
values = apal[c(4,3,2,1)]) +
scale_y_continuous(breaks = c(0, 20, 40, 60, 80, 100), labels = paste0(c(0, 20, 40, 60, 80, 100), '%')) +
labs(x = 'Variance (denominator)', y = '', fill = 'Variance source') +
jftheme
}
plot_predictions <- function(amodel, point_alpha){
adt <- as.data.table(amodel@frame)
adt[, wave := factor(as.numeric(gsub('ses-(\\d+)', '\\1', sess)))]
newdata <- unique(adt, by = c('id', 'sess', 'wave'))[, c('H', 'z') := list(0, NULL)]
newdata$z_prime <- predict(amodel, newdata = newdata)
newdata$z_prime_id <- predict(amodel, newdata = newdata, re.form = ~(1|id))
ggplot(newdata, aes(x = wave, y = tanh(z_prime), group = id)) +
#geom_violin(data = adt, aes(group = NULL, y = tanh(z)), size = 0, alpha = .5, fill = apal[[1]]) +
geom_point(data = adt, aes(group = NULL, y = tanh(z)), alpha = point_alpha, color = apal[[1]], position = position_jitter(), size = .5) +
geom_point(color = apal[[3]], alpha = 1) +
geom_line(color = apal[[3]]) +
geom_rug(data = unique(newdata, by = 'id'), aes(y = z_prime_id, x = NULL), color = apal[[4]], alpha = 1, sides = 'tr', length = unit(0.03, "npc")) +
geom_hline(yintercept = 0, color = apal[[1]]) +
coord_cartesian(ylim = c(-.025, .65), xlim = c(1, 10)) +
labs(y = 'FC correlation', x = 'Session') +
jftheme
}
max_points <- unlist(sea_schaefer400_7_withinnet[, .N, by = net1][, list(max = max(N))])
point_alphas <- sea_schaefer400_7_withinnet[, .N, by = net1][, p := 1 - N / max_points][, pomp := (p - min(p)) / (max(p) - min(p))][, alpha := .025 + pomp*.1]
font_add_google("Didact Gothic", "Didact Gothic")
showtext_auto()
for(i in 1:length(model_fits)){
model_fit <- model_fits[[i]]
network_name <- networks[[i]]
cat('\n\n### ', network_name, '{.tabset}\n\n')
cat('\n\n#### Plot\n\n')
point_alpha <- point_alphas[net1 == network_name, alpha]
print(plot_predictions(model_fit, point_alpha = point_alpha) + plot_percents(proportion_variance_tables[[i]]) +
plot_layout(ncol = 2, widths = c(4,1)))
cat('\n\n#### Table (%)\n\n')
print(knitr::kable(tidyr::spread(proportion_variance_tables[[i]], out_of, percent), digits = 1))
cat('\n\n#### Model Summary\n\n')
cat('\n```\n')
print(summary(model_fit))
cat('\n```\n')
}Cont
Plot
Table (%)
| source | total | RE |
|---|---|---|
| ID | 0.0 | 0.1 |
| ID/Session | 5.1 | 99.9 |
| Residual | 94.9 | NA |
Model Summary
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 43287.4
Scaled residuals:
Min 1Q Median 3Q Max
-5.1523 -0.6818 -0.0690 0.6212 6.1048
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 3.936e-03 0.062736
id (Intercept) 5.080e-06 0.002254
Residual 7.345e-02 0.271023
Number of obs: 186298, groups: sess:id, 289; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.219505 0.003766 58.28
Default
Plot
Table (%)
| source | total | RE |
|---|---|---|
| ID | 0.0 | 0 |
| ID/Session | 3.8 | 100 |
| Residual | 96.2 | NA |
Model Summary
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 186542
Scaled residuals:
Min 1Q Median 3Q Max
-4.8942 -0.6869 -0.0774 0.6111 6.0356
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.003199 0.05656
id (Intercept) 0.000000 0.00000
Residual 0.080268 0.28332
Number of obs: 587200, groups: sess:id, 289; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.219510 0.003348 65.57
convergence code: 0
boundary (singular) fit: see ?isSingular
DorsAttn
Plot
Table (%)
| source | total | RE |
|---|---|---|
| ID | 0.3 | 3.8 |
| ID/Session | 6.5 | 96.2 |
| Residual | 93.3 | NA |
Model Summary
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 58752.1
Scaled residuals:
Min 1Q Median 3Q Max
-5.3612 -0.6951 -0.0866 0.6100 5.6393
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.0061258 0.07827
id (Intercept) 0.0002388 0.01545
Residual 0.0879342 0.29654
Number of obs: 141896, groups: sess:id, 289; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.274858 0.005462 50.32
Limbic
Plot
Table (%)
| source | total | RE |
|---|---|---|
| ID | 0.3 | 3.1 |
| ID/Session | 9.9 | 96.9 |
| Residual | 89.8 | NA |
Model Summary
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: -4339.2
Scaled residuals:
Min 1Q Median 3Q Max
-4.4720 -0.6583 -0.0461 0.6085 5.7245
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.0057390 0.07576
id (Intercept) 0.0001826 0.01351
Residual 0.0519942 0.22802
Number of obs: 43652, groups: sess:id, 289; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.195606 0.005212 37.53
SalVentAttn
Plot
Table (%)
| source | total | RE |
|---|---|---|
| ID | 0.0 | 0 |
| ID/Session | 6.2 | 100 |
| Residual | 93.8 | NA |
Model Summary
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 24328.5
Scaled residuals:
Min 1Q Median 3Q Max
-6.1158 -0.6666 -0.0501 0.6034 6.4400
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.004479 0.06693
id (Intercept) 0.000000 0.00000
Residual 0.068308 0.26136
Number of obs: 151093, groups: sess:id, 289; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.262361 0.003994 65.69
convergence code: 0
boundary (singular) fit: see ?isSingular
SomMot
Plot
Table (%)
| source | total | RE |
|---|---|---|
| ID | 0 | 0 |
| ID/Session | 10 | 100 |
| Residual | 90 | NA |
Model Summary
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 98364
Scaled residuals:
Min 1Q Median 3Q Max
-7.0458 -0.6704 -0.0996 0.5656 7.3660
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 0.008267 0.09092
id (Intercept) 0.000000 0.00000
Residual 0.074277 0.27254
Number of obs: 407269, groups: sess:id, 289; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.305983 0.005366 57.02
convergence code: 0
boundary (singular) fit: see ?isSingular
Vis
Plot
Table (%)
| source | total | RE |
|---|---|---|
| ID | 0.0 | 0.8 |
| ID/Session | 5.2 | 99.2 |
| Residual | 94.8 | NA |
Model Summary
Linear mixed model fit by REML ['lmerMod']
Formula: z ~ 1 + (1 | id/sess)
Data: this_net_dt
REML criterion at convergence: 168016.6
Scaled residuals:
Min 1Q Median 3Q Max
-3.5945 -0.7092 -0.0977 0.6328 4.8807
Random effects:
Groups Name Variance Std.Dev.
sess:id (Intercept) 6.155e-03 0.078456
id (Intercept) 4.922e-05 0.007015
Residual 1.128e-01 0.335917
Number of obs: 254358, groups: sess:id, 289; id, 30
Fixed effects:
Estimate Std. Error t value
(Intercept) 0.317729 0.004837 65.69